#Load data
GetData<-function(size=5000,boot=TRUE,directory="",y0=1986,y1=2015){
  if(directory==""){
    Datawd<-getwd()
  }else{
    Datawd<-directory
  }
  yy0<-substr(y0,(nchar(y0)+1)-2,nchar(y0))
  yy1<-substr(y1,(nchar(y1)+1)-2,nchar(y1))
  #Load data for years y0 and y1
  load(file.path(Datawd,paste0("Data",yy0,"-",yy1,".Rda")))
  
  vars<-c("lrwage3","union","pubsect","manuf","nonwhite",
          "female","partte","married","smsa",
          "exper" , "exper2" ,
          "exp_cl" , 
          paste0("expe",1:length(unique(Data$exp_cl))), 
          "years_educ" , "educ_cl" , 
          paste0("educ",1:length(unique(Data$educ_cl))), 
           "edex" ,
          "region",
          paste0("reg",1:length(unique(Data$region))),
          "ee_cl",
          paste0("ee",1:length(unique(Data$ee_cl))),
          "ind21",
          paste0("indu",1:length(unique(Data$ind21))),
          "state",
          paste0("state",1:length(unique(Data$state))),
          "orgwgt","year")
  
  Z0<-Data[Data$year==y0,vars]
  Z1<-Data[Data$year==y1,vars]
  if (boot==TRUE){
    w="orgwgt"
    probab0=eval(parse(text=paste("Z0$",w,"/sum(Z0$",w,")",sep="")))
    set.seed(1969)
    i0 <- sample(x=1:nrow(Z0),size=size, replace = FALSE,prob=probab0)
    Z0<-Z0[i0,vars]
        
    probab1=eval(parse(text=paste("Z1$",w,"/sum(Z1$",w,")",sep="")))
    set.seed(1983)
    i1 <- sample(x=1:nrow(Z1),size=size, replace = FALSE,prob=probab1)
    Z1<-Z1[i1,vars]
  }
  
  rm(Data)
  
  return(list(Z0=Z0,Z1=Z1))
}

#Function to get coefficients from the QR
ParRq<-function(t,f=formula,d=Z0,w="orgwgt",mth="pfn"){
  require(quantreg)
  if (w!=""){
    h<-eval(parse(text=paste("rq(f, data=d,tau=t,weights=",w,",method='",mth,"')$coef",sep="")))
  } else{
    h<-rq(f, data=d,tau=t,method=mth)$coef
  }
  return(h)
}

#Function to get the coefficients from the QR in parallel
ParCoefRq<-function(m=4500,formula,d,w="orgwgt",mth="pfn",seed=13,taus=NULL,nodes=15,ParRq){
  require(doParallel)
  cl <- makeCluster(nodes, type = "SOCK")
  registerDoParallel(cl)
  set.seed(seed)
  if (is.null(taus)){
    eps<-1/70
    u<-runif(m, min = 0+eps, max = 1-eps)
    CoefRq<-foreach(tt=u,.combine='rbind',.inorder=FALSE)%dopar%{
      ParRq(tt,f=formula,d=d,w=w,mth=mth)
    }
  }else{
    CoefRq<-foreach(tt=taus,.combine='cbind',.inorder=TRUE)%dopar%{
      ParRq(tt,f=formula,d=d,w=w,mth=mth) 
    }
  }
  stopCluster(cl)
  return(CoefRq)
}

SampW<-function(Z,CoefRq,w="",m=4500,seed=14){
  set.seed(seed)
  if(w==""){
    ii <- sample(x=1:nrow(Z),size=m, replace = TRUE)
  }else{
    probab=eval(parse(text=paste("Z$",w,"/sum(Z$",w,")",sep="")))
    ii <- sample(x=1:nrow(Z),size=m, replace = TRUE,prob=probab)
  }
  Ze<-Z[ii,]
  rownames(Ze)<-seq(1,m,by=1)
  Ze<-Ze[,colnames(Ze)%in%colnames(CoefRq)]
  Ze<-cbind(rep(1,m),Ze)
  names(Ze)[1]<-colnames(CoefRq)[1]
  return(list(We=rowSums(Ze*CoefRq),Ze=Ze))
}

SampWcf<-function(SW1e,I1,I0,seeds,m,ss=T){
  #I1 is the index of classes C(1),..,C(J) in sample SW1e$Ze at t=1
  #I0 is the index of classes C(1),..,C(J) in populaton Z0 at t=0
  #this should be true: 
  #length(I1)==length(I0)
  #m is the sample size of SW1e$Ze
  #for reproducibility set seeds
  if(length(I1)!=length(I0)){
    stop("number of classes in t=1 and t=0 is different")
  }else{
    temp<-NULL
    for (c in 1:length(I1)){
      #Subsample of W1e's of class c storaged in WIc
      t<-paste("WI",c,"<-SW1e$We[I1$i",c,"]",sep="")
      eval(parse(text=t))
      t<-paste("names(WI",c,")<-1:length(WI",c,")",sep="")
      eval(parse(text=t))
      
      #Frequency of class c in the original data at t=0
      t<-paste("fi0<-sum(I0$i",c,",na.rm=T)/length(I0$i",c,")",sep="")
      eval(parse(text=t))
      #Sample size to generate
      mm<-round(m*fi0)
      if (ss==FALSE){
        t<-paste("ii<-sample(x=1:length(WI",c,"),size=mm, replace = TRUE)",sep="")
        eval(parse(text=t))
      }else{
        set.seed(seeds[c])
        t<-paste("ii<-sample(x=1:length(WI",c,"),size=mm, replace = TRUE)",sep="")
        eval(parse(text=t))
      }
      t<-paste("temp<-c(temp,WI",c,"[ii])",sep="")
      eval(parse(text=t))
    }
    names(temp)<-1:length(temp)
    return(temp)
  } 
}

BootW<-function(W,seed=69){
  set.seed(seed)
  ii<-sample(x=1:length(W),size=length(W), replace = TRUE)
  Bw<-W[ii]
  return(Bw)
}

change<-function(W0,W1,probs=c(0.01,0.1,0.25,0.5,0.75,0.9,0.99),seed=69,Boot=TRUE){
  require(reldist)
  if(Boot==TRUE){
    w0<-BootW(W0,seed=seed)
    w0<-w0[!is.na(w0)]
    w1<-BootW(W1,seed=seed)
    w1<-w1[!is.na(w1)]
  }else{
    w0<-W0[!is.na(W0)]
    w1<-W1[!is.na(W1)]
  }
  
  M0<-quantile(w0,probs=probs,na.rm=TRUE)
  M1<-quantile(w1,probs=probs,na.rm=TRUE)
  
  M0<-c(M0,gini(w0))
  M1<-c(M1,gini(w1))
  names(M0)[length(names(M0))]="GiniLnW"
  names(M1)[length(names(M1))]="GiniLnW"
  
  M0<-c(M0,gini(exp(w0)))
  M1<-c(M1,gini(exp(w1)))
  names(M0)[length(names(M0))]="GiniW"
  names(M1)[length(names(M1))]="GiniW"
  
  return(M1-M0)
}

#Function get bootstrap confidence interval in parallel
ParBootCI<-function(R=1000,nodes=15,BootW,fun=change,...){
  require(doParallel)
  cl <- makeCluster(nodes, type = "SOCK")
  registerDoParallel(cl)
  
  Boot<-foreach(r=1:R,.combine='rbind',.inorder=FALSE)%dopar%{
    BootW<-BootW
    fun(...,seed=r)
  }
  stopCluster(cl)
  return(Boot)
}


actual<-function(W,probs=c(0.01,0.1,0.25,0.5,0.75,0.9,0.99),seed=69,Boot=TRUE){
  require(reldist)
  if(Boot==TRUE){
    w<-BootW(W,seed=seed)
  }else{
    w<-W
  }
  M<-quantile(w,probs=probs,na.rm=TRUE)
  M<-c(M,gini(w))
  names(M)[length(names(M))]="GiniLnW"
  M<-c(M,gini(exp(w)))
  names(M)[length(names(M))]="GiniW"
  return(M)
}

Ind_coef<-function(X0,X1,coef){
  temp<-X1
  temp[,coef]=X0[,coef]
  return(temp)
}

Tab<-function(T1,T2,T3){
  T1[rownames(T1)%in%c("GiniW","GiniLnW"),]<-100*T1[rownames(T1)%in%c("GiniW","GiniLnW"),]
  T2[rownames(T1)%in%c("GiniW","GiniLnW"),]<-100*T2[rownames(T1)%in%c("GiniW","GiniLnW"),]
  
  Table<-matrix(rep("",3*nrow(T1)*ncol(T1)),nrow=3*nrow(T1),ncol=ncol(T1))
  colnames(Table)<-colnames(T1)
  rnames<-matrix(rep("",3*nrow(T1)),nrow=3*nrow(T1),ncol=1)
  rnames[seq(1,3*nrow(T1),by=3)]<-rownames(T1)
  rownames(Table)<-rnames
  
  Table[seq(1,3*nrow(T1),by=3),]<-format(round(T1,digits=3),scientific=FALSE)
  CIL<-matrix(rep("",nrow(T2)*ncol(T2)*0.5),nrow=nrow(T2),ncol=ncol(T2)*0.5)
  colnames(CIL)<-colnames(Table)[3:length(colnames(Table))]
  rownames(CIL)<-rownames(T2)
  
  low<-T2[,colnames(T2)%in%"2.5%"]
  ci_ni<-T2[!rownames(T2)%in%c("GiniW","GiniLnW"),colnames(T2)%in%"2.5%"]<=0
  CIL[rownames(CIL)%in%rownames(ci_ni),][ci_ni]<-paste(format(round(T2[rownames(T2)%in%rownames(ci_ni),colnames(T2)%in%"2.5%"][ci_ni],digits=3),scientific=FALSE),";",sep="")
  CIL[rownames(CIL)%in%rownames(ci_ni),][!ci_ni]<-paste(format(round(T2[rownames(T2)%in%rownames(ci_ni),colnames(T2)%in%"2.5%"][!ci_ni],digits=3),scientific=FALSE),";",sep="")
  
  ci_ni<-T2[rownames(T2)%in%c("GiniW","GiniLnW"),colnames(T2)%in%"2.5%"]<0
  CIL[rownames(CIL)%in%rownames(ci_ni),][ci_ni]<-paste(format(round(T2[rownames(T2)%in%rownames(ci_ni),colnames(T2)%in%"2.5%"][ci_ni],digits=3),scientific=FALSE),";",sep="")
  CIL[rownames(CIL)%in%rownames(ci_ni),][!ci_ni]<-paste(format(round(T2[rownames(T2)%in%rownames(ci_ni),colnames(T2)%in%"2.5%"][!ci_ni],digits=3),scientific=FALSE),";",sep="")
  
  CIU<-matrix(rep("",nrow(T2)*ncol(T2)*0.5),nrow=nrow(T2),ncol=ncol(T2)*0.5)
  colnames(CIU)<-colnames(Table)[3:length(colnames(Table))]
  rownames(CIU)<-rownames(T2)
  
  high<-T2[,colnames(T2)%in%"97.5%"]
  ci_ni<-T2[!rownames(T2)%in%c("GiniW","GiniLnW"),colnames(T2)%in%"97.5%"]<=0
  CIU[rownames(CIU)%in%rownames(ci_ni),][ci_ni]<-format(round(T2[rownames(T2)%in%rownames(ci_ni),colnames(T2)%in%"97.5%"][ci_ni],digits=3),scientific=FALSE)
  CIU[rownames(CIU)%in%rownames(ci_ni),][!ci_ni]<-format(round(T2[rownames(T2)%in%rownames(ci_ni),colnames(T2)%in%"97.5%"][!ci_ni],digits=3),scientific=FALSE)
  
  ci_ni<-T2[rownames(T2)%in%c("GiniW","GiniLnW"),colnames(T2)%in%"97.5%"]<=0
  CIU[rownames(CIU)%in%rownames(ci_ni),][ci_ni]<-format(round(T2[rownames(T2)%in%rownames(ci_ni),colnames(T2)%in%"97.5%"][ci_ni],digits=3),scientific=FALSE)
  CIU[rownames(CIU)%in%rownames(ci_ni),][!ci_ni]<-format(round(T2[rownames(T2)%in%rownames(ci_ni),colnames(T2)%in%"97.5%"][!ci_ni],digits=3),scientific=FALSE)
  
  CI<-matrix(paste(CIL,CIU,sep=""),nrow=nrow(T2),ncol=ncol(T2)*0.5)
  rownames(CI)<-rownames(CIL)
  colnames(CI)<-colnames(CIL)
  
  star<-matrix(rep("",nrow(low)*ncol(low)) ,nrow=nrow(low),ncol=ncol(low))
  star[!(low<0 & high >0)]<-"*"
  colnames(star)<-colnames(CI)
  rownames(star)<-rownames(CI)
  
  Table[seq(1,3*nrow(T1),by=3),colnames(Table)%in%colnames(star)]<-paste(Table[seq(1,3*nrow(T1),by=3),colnames(Table)%in%colnames(star)],star,sep="")
  Table[seq(2,3*nrow(T1),by=3),colnames(Table)%in%colnames(CI)]<-CI
  Table[seq(3,3*nrow(T1),by=3),4:length(colnames(Table))]<-format(round(T3,digits=3),scientific=FALSE)
  #format(Table,justify="centre")
  
  Res<-matrix(rep("",3*nrow(T1)),nrow=3*nrow(T1),ncol=1)
  Res[seq(1,3*nrow(T1),by=3),]<-format(round(T1[,"Marg"]-T1[,"Cova"]-T1[,"Coef"],4),scientific=FALSE)
  Res[seq(3,3*nrow(T1),by=3),]<-format(round((T1[,"Marg"]-T1[,"Cova"]-T1[,"Coef"])/T1[,"Marg"],3),scientific=FALSE)
  
  Table<-cbind(Table[,1:5],Res,Table[,6:ncol(Table)])
  colnames(Table)[colnames(Table)==""]<-"Res"
  return(Table)
}

save.xlsx <- function (file, ...){
  require(xlsx, quietly = TRUE)
  objects <- list(...)
  fargs <- as.list(match.call(expand.dots = TRUE))
  objnames <- as.character(fargs)[-c(1, 2)]
  nobjects <- length(objects)
  for (i in 1:nobjects) {
    if (i == 1)
      write.xlsx(objects[[i]], file, sheetName = objnames[i])
    else write.xlsx(objects[[i]], file, sheetName = objnames[i],
                    append = TRUE)
  }
  print(paste("Workbook", file, "has", nobjects, "worksheets."))
}

#Function to get time attribute
`%.%` <- function(o, a) attr(o, as.character(substitute(a)))
